Random Networks of Spiking Neurons: 
Instability in the Xenopus tadpole moto-neural pattern 
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A large network of integrate-and-fire (IF) neurons is studied analytically when the synaptic 
weights are independently randomly distributed according to a gaussian distribution with arbitrary 
mean. The relevant order parameters are identified, and it is shown that such network is statistically 
equivalent to an ensemble of independent IF neurons with each input signal given by the sum of a 
self-interaction deterministic term and a gaussian coloured noise. The model is able to reproduce 
the quasi-synchronous oscillations, and the dropout of their frequency, of the CNS neurons of the 
swimming Xenopus tadpole. The cause of the instability of low firing-rate self-sustained activity is 
identified and explained. Predictions from the model are proposed for future experiments. 
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In the recent past, the study of neuronal dynamics has 
been a major topic, both for the importance it has in neu- 
roscience and for the challenge it represents to modern 
non-linear mathematics. Particular interest is evoked by 
the recurrent networks of many neurons, since they are 
present in many neuronal systems, in mammal cortex 
as well as in small vertebrates' central nervous system. 
Specific of the latter is the case of the Xenopus tadpole 
embryo, in which the relationship between sensorial stim- 
uli and the consequent behaviour of the animal can be 
studied in detail at the neuronal level 0J^. Pools of recip- 
rocally connected excitatory neurons, distributed bilater- 
ally along the spinal cord, are the elements of the moto- 
neural circuits of the Xenopus tadpole that are mainly 
responsible for the generation of the swimming pattern, 
consisting of regular oscillations of the body of the tad- 
pole on the horizontal plane. Each pool includes a few 
hundreds of excitatory neurons receiving inputs from the 
sensory system, and producing outputs that are sent to 
the other neurons of the same pool, to pools of inhibitory 
neurons that act on the symmetrical contra-lateral exci- 
tatory pool, and to neurons responsible for muscular ac- 
tivation The only inhibitory signals to the neurons 
of the pool come from outside, mainly caused by activ- 
ity of the contra-lateral pool. The rhythm of the swim 
may be not, or at least not only, due to the reciprocal 
inhibition between the pools on the opposite sides of the 
spinal cord. Indeed, it has been shown that any exci- 
tatory pool is able to present a series of oscillations even 
if the commissural connections between the two sides of 
the spinal system are removed and the inhibition is sup- 
pressed. The mechanism producing the oscillations in 
this isolated pool is not fully understood yet. 

In order to investigate on the basic mechanisms of the 
neuronal activity in recurrent networks, particularly in 
the case of the Xenopus, in the present work a large net- 
work of N integrate-and-fire (IF) neurons with gaussian 
synapses is studied. The membrane potentials are here 



referred to the reset value = 0, and are measured in 
terms of the threshold potential Vth — 1- Time is mea- 
sured in units of membrane time constant r = 1. The 
equation followed by the membrane potential of neuron 
i between any two consecutive spikes is 
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dt 
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where / is the external current, and T™ is the instant 
in which neuron j fires its m-th spike. The coefficient 
Wij is the synaptic weight that multiplies the ionic cur- 
rent density due to the spikes arriving to neuron i from 
neuron j. The weights {wij,i ^ j} are taken to be in- 
dependently randomly distributed according to a gaus- 
sian p.d.f. with mean and variance equal to /x/iV and 
CT^/iV, respectively. Thus, every Wij includes the 'sign' 
of the interaction: Positive for an excitatory synapse, 
and negative for an inhibitory one. The weights are 
quenched and are not required to be symmetric. The 
function 3(t) contains the details of the ionic current den- 
sity through the synaptic channels, and in all the analyt- 
ical part of the work its functional form is not relevant. 
In the numerical study, it is taken as a a-function that 
accounts for finite time dynamics of the synaptic chan- 
nels: J{t) = o?{t — Ta)9{t — Ta) exp{—a{t — To)}, where 
also the axonal delay time Tq, for simplicity taken to be 
the same for all the axonal projections, has been included 
to simplify the notations. 

Let T be the matrix of the firing times of all the neu- 
rons. Any quantity relative to the network must be a 
function of T, for example F{T). Of particular interest 
are the quantities that do not depend on the specific re- 
alization of the network when the number of neurons is 
very large (thermodynamic limit). This property (self- 
averaging) allows one to substitute the calculation of F 
for a specific realization, usually unfeasible, with the cal- 
culation of its average over all the possible realizations 
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of the random network. Normally, one cannot prove a 
priori if the considered quantity self-averages, and has to 
calculate also the fluctuation of F across the realizations. 
The average of F across the ensemble of networks is 

< F{T) >= j dw V{w) j dT F{T)^6{G{w,T)), 

(2) 

where the partition function Z = J dT_ 5 {G{w, T)) guar- 



antees the correct normalization of the (5-function. The 
deterministic dynamics of the system is represented by 
the matrix G(w, T), whose generic element gives a re- 
lation between two consecutive spikes of a neuron and 
the spikes previously fired by the other neurons (coming 
from the IF dynamics, like in |^). For the application of 
interest in this letter, the uniform initial conditions are 
adopted here according to which all the neurons are ini- 
tially stably at membrane rest potential. Thus, for any 
n and i one has: 
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dt 



J2jit-Tp). 



(3) 



With simple appropriate modifications, the mathemati- 
cal framework here developed can be applied also to many 
cases in which the neurons of the net are given different 
initial conditions and different external currents. 

The quantities of interest in this paper can be obtained 
through the usual techniques based on derivatives of ficti- 
tious fields introduced at a later stage during the calcula- 
tion of the 'free-energy' ^ In Z 3>, and thus the following 
analysis concentrates on that. 

To calculate the free-energy, the replica trick is used: 

< In Z hm - In < Z' >, (4) 

I — tO r 

where the limit is obtained after having calculated the 
average for integer r and then having prolonged <C 
Z^ ^ analytically to real exponents. Substituting the 
(5-functions with their complex Fourier expansions, one 
has that 

<Z'->= J DwDTDs exp{ ^ is"'"Gr(w,T")}, 

(5) 

being z, n, a indices of, respectively, neuron, spike, and 
replica. Through functional integration over one-time 
and two-times auxiliary functions one can decouple neu- 
rons. Due to the absence of dynamic noise and of any 
form of average over the initial conditions, the replicas 
are identical also in their kinematics, and this allows for 
an exact replica-symmetry assumption. Following the 
symmetry, after the saddle-point evaluation allowed by 
the extensiveness of the network {N —^ oo), only two 
order parameters remain in the formulas, both having 
direct physical meaning as averages over the ensemble of 
networks. They are: 

i m 

and 



i m n 

(7) 

Thus, the order parameters 'filter' the spike train 
{r",Vn} fired by neuron i through the current density 
function J(-). It suggests that the correct smoothing 
function to analyse the information conveyed by a spike 
train in real data may be the synaptic function J . It can 
be shown that the arithmetic averages over the neurons, 
inside the double angle brackets, have vanishing fluctua- 
tions in the thermodynamic limit of the mean-field theory 
(self-averaging); then the order parameter physics make 
sense also for any single network out of the ensemble. At 
the level of the single network, the parameter v{t) may 
be seen as an estimate of the average firing rate in the 
network if the firing rate is measured in a delayed back- 
time window about 3/a long. Actually, v{t) weights the 
preceding spikes according to their synaptic 'efficacy' at 
time t. The order parameter q{ti,t2) is the average auto- 
correlation of the 'filtered' spike trains from the neurons 
of the network. Another two-times physical order pa- 
rameter would appear if dynamic noise was present in 
the network. In fact, both their physical meaning and 
the analysis of the saddle-point equations show that the 
two two-times order parameters are equal to each other 
thanks to the absence of dynamic noise. 

The analysis allows for the interpretation of the mean- 
field system in terms of an independent IF neuron with 
total input equal to I+T{t)+^ ^(i)j being f(i) a gaussian 
stochastic process with auto-correlation a^q{t,t'). Indi- 
cating the average over the realizations of the stochastic 
input with single angle brackets, the order parameters 
follow the equations 

:^it)=<Y,J{t~Tn>, (8) 



2 



and 

q{t, t') =< Jit - T") ^ J{t' ~ T") >, (9) 

m n 

where T" is the n-th spike time of the equivalent neu- 
ron, and the sequence of spiking times is determined 
by the IF dynamics of the equivalent neuron with input 
I + f{t)+iiiy{t). 

Solving analytically this new set of equations, as well 
as the original saddle-point equations, is a difficult task. 
However, simulating the equivalent neuron is easy nu- 
merically, and finite-size effects are no longer present, 
as noted by |^ and since the thermodynamic limit 
has already been performed. In order to simulate the 
equivalent system and obtain quantities of interest for 
the original network, use is made here of the technique 
introduced by j^, where the correlated stochastic input 
is constructed during the running of a quite large number 
of independent simulations of the equivalent neuron. In 
this method, the correlation function is estimated from 
sample averages, and from that the gaussian coloured 
signal is created. As suggested by since each one- 

dimensional simulation is independent from the others, 
the error on the averages across a set of M trials should 
scale as I/a/M- 

An exhaustive analysis in the parameters' space is be- 
ing worked out, and will be presented elsewhere, as well 
as an analysis of the possibility of substituting the aver- 
ages over the realizations of the gaussian noise in eqs. ^ 
and H with time averages over a single realization of the 
noise (ergodicity) that may be useful, e.g., in numerical 
computations. The rest of this letter concentrates on the 
particular range of parameters that is of interest for the 
study of the neural control of the swimming pattern in 
the simple vertebrate Xenopus tadpole. 

It has been shown |^ that any excitatory pool surgi- 
cally isolated from the rest of the tadpole nervous system 
is able to perform several oscillation cycles after a brief 
external stimulation that is provided to a large number 
of neurons inside the pool. It is biologically plausible to 
assume that the synapses between neurons of the excita- 
tory pool in the Xenopus spinal system have weights very 
similar to each other (in the model it means: /i > and 
cr ^ /i). In the experiments on the tadpole considered 
here a brief external input is provided simultaneously to 
a large number of neurons in a quiescent pool. Thus, the 
neurons start firing synchronously. The external stim- 
ulus is brief, and after its removal the network cannot 
stay stable. Simulations of the equivalent neuron show 
that the already present activity decays to zero if fi is 
not too large; if /i is above a certain bound, the activ- 
ity in the network increases restlessly to saturation. The 
reason of such instability is discussed later in this let- 
ter. For the moment, consider a case with /i below the 
aforementioned critical bound. Figure ^ shows the time 
evolution of the order parameter i^it) in a network with 
fj, = .9 and — .001. At the beginning, all the neurons 
are in a quiescent state, and a constant uniform input 



current / is applied (0*'' time-step). After a while, the 
external source is abruptly removed, and the network is 
let to its own dynamic evolution, during which it exhibits 
a quasi-synchronous activity whilst the maximal magni- 
tude of the order parameter i/(i) decays, correspondingly 
to the also evident decrease in spike frequency. Evident 
as well is the progressive loss of synchrony (vertical bars; 
cf. figure caption). The frequency drop continues until 
the network rests in a quiescent state. Here a case of 
quite rapid dropout is shown, but larger /x (still below 
the bound) can provide slower dropout. The smooth- 
ness of the curve near its local minima is an effect of the 
small variance of the synaptic strengths, that perturbs 
the synchrony of the spikes, with a consequent smooth- 
ing of the average. In complete absence of noise, that 
is when all the synapses have equal weight, the curve 
of 1/(1) presents sharp downward peaks in the local min- 
ima (discontinuities of the first time derivative) in de- 
layed correspondence with the spike events. The decrease 
in synchrony may account for progressive single-neuron 
spike dropout found in experiments, often interpreted as 
a spike or synaptic failure. Increasing a causes quicker 
synchrony loss. These results reproduce the experimen- 
tal data of frequency dropout Q]. A more quantitative 
comparison will be presented when appropriate experi- 
mental data are available. Meanwhile, direct simulations 
of the IF network have been carried out: the agreement 
between theoretical predictions and simulation outputs 
is very good already for networks with as few as 100 neu- 
rons. This also supports the use of the thermodynamic 
limit as an approximation of the real finite network. 

The form of instability of self-sustained activity pre- 
sented by the network in the region of the parameters' 
space that mimics the tadpole biological values, is due to 
the absence of any synaptic currents' and channels' reset 
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t (in time-steps) 

FIG. 1. Time evolution of the order parameter u{t) in the net- 
work that mimics a Xenopus tadpole excitatory pool of motor neu- 
rons. The first three arrows indicate the spikes due to the initial 
external input. The fourth arrow (below the abscissae axis) indi- 
cates the instant of the removal of the external source. The bars 
indicate the standard deviations of the average spike trains (that is, 
±y^q{t, t) — v(t)v(t)); they clearly report a progressive synchrony 
loss, {^l = .9, 0-2 = 001, a = 5., Ta = .2, I = 15 (only until the 
25*'' time-step), M = 10000, time-step=.01.) 
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after the generation of a spike, so that they affect the 
dynamics of the membrane independently of it. It is eas- 
ily understood how this causes instability in the isolated 
network: Consider the case of null synaptic variance, for 
simplicity, and no external input, and assume that all 
the neurons fire simultaneously the first spike, as an ini- 
tial condition, for example, so that there are no resid- 
ual synaptic currents from previous spikes or external 
sources. Thus, after time Tq, any neuron of the network 
receives N spikes, and starts integrating the synaptic cur- 
rents. With some analysis of the equation it is possible 
to show that the minimum synaptic strength fi that al- 
lows the membrane potential to reach the threshold de- 
pends on a as depicted in figure |[ If /i is below the 
curve, the neurons do not fire any other spike. On the 
contrary, if fi is above the curve the neurons fire at finite 
time from the arrival of the first spikes; after time r^, 
every neuron receives the second wave of spikes, and the 
new synaptic currents are integrated together with the 
remains of the currents generated by the first wave, and 
the neuron fires again. This process is iterated, and the 
effect of the remaining currents is to decrease monoton- 
ically the inter-spike interval, thus driving the network 
to its maximum firing rate. Simulations of the equiva- 
lent neuron seem to indicate that the introduction of a 
small variance in the synaptic weights lowers the critical 
bound. 

To test the hypothesis about the network instability, 
numerical solutions have been found introducing in the 
model the reset of the input currents to zero immedi- 
ately after every generated spike. The calculations as 
presented here cannot include this effect, but for small 
values of cr^, like in the case of the tadpole modelling, the 
results are still valid for not too long times. Simulations 
with this proviso show that when ^ is above the bound, 
the synaptic variance only destroys synchrony, while the 
network is in a stable firing state. When /i is below the 
critical bound, the network is not able to fire more than 
one spike after the removal of the external stimulus. One 
can also obtain stability at low firing-rates substituting 
the synaptic a-function with a function that vanishes at 
finite time, without a currents' hard reset. If two con- 
secutive spikes do not come at intervals shorter than the 
time necessary for the synaptic current to vanish, current 
remains do not cumulate and the network can fire stably 
for large ^. This latter case would be easily and exactly 
reproducible in the mathematical framework developed 
here. 

Since the synaptic currents in the Xenopus tadpole em- 
bryo seem to be well approximated by the a-function , 
the results of the present work predict that the mean 
synaptic strength ^ in the real network is below the crit- 
ical bound, and the series of oscillations subsequent to 



the removal of the external stimulus is due to a 'surfing' 
on the remains of synaptic currents until their exhaus- 
tion. The present model also predicts a progressive loss 
of synchrony during the oscillatory transient, due to the 
small deviation of the synaptic strengths from the mean. 
A further implication of the present results is that some 
mechanisms proposed in the past as responsible for the 
oscillation dropout, like synaptic depression or habitua- 
tion, are not needed. Besides, the present work also gives 
a proof of the instability of the synchronous state in iso- 
lated fully-connected IF random neural networks with 
< cr < /i. 

If the predictions about mean synaptic strength, syn- 
chrony loss, and single-spike dropout are verified experi- 
mentally, the present model may be useful in understand- 
ing some neural mechanisms subserving the small varte- 
brate's motion. With the aid of the mathematical results 
presented here, other, more general topics are under cur- 
rent study, like the possible presence of several attrac- 
tors, chaos, and existence of a glassy phase in random 
networks. 

I thank Prof. Paul Bressloff for having brought the 
spiking neurons' random networks to my attention, and 
for useful discussions. 
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FIG. 2. Minimum synaptic strength necessary to generate the 
second spike (cf. text), as a function of the synaptic constant a. 
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